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Abstract 

To avoid the combinatorial computational cost of configuration interaction (Cl), we have pre¬ 
viously introduced the symmetric tensor decomposition Cl (STD-CI) method, where we take ad¬ 
vantage of the antisymmetric nature of the electronic wave function and express the Cl coefficients 
compactly as a series of Kronecker product states (STD series) [W. Uemura and O. Sugino, Phys. 
Rev. Lett. 109, 253001 (2012)]. Here we extend the variational degrees of freedom by using dif¬ 
ferent molecular orbitals for different terms in the STD series. This scheme is equivalent to the 
linear combination of the Hartree-Fock-Bogoliubov state or the antisymmetrized geminal powers 
(AGP). The total energy converges very rapidly within 0.72 /xHartree taking only 10 terms for the 
water molecule, and the convergence is likewise fast for Hubbard tetramers. The computational 
cost scales as the fifth power of the number of electrons and the square of the number of terms in 
the STD series, indicating the promise of this AGP-based scheme for highly accurate and efficient 
computation of quantum systems. 
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I. INTRODUCTION 


Determining the ground state of a many-body system is the most basic problem in many 
helds of science, such as condensed matter physics, nuclear physics and quantum chem¬ 
istry. Many numerical approaches have been developed so far to obtain the wavefunction. 


In the conhguration interaction (CI)[l|, multiple quasiparticle states are generated from a 
mean-held wavefunction and the coefficients for their linear combination are determined by 
solving an eigenvalue problem, while in the variational Monte Carlo (VMC) the quasi¬ 
particles are generated stochastically to obtain the expectation value of the total energy 
to be variationally determined. In the tensor network (TN) framework[^, the total energy 
is variationally determined using specihc quasiparticle states generated according to an as¬ 
sumed TN. In those popular schemes, accessible degrees of freedom are usually not very 
large. This is particularly the case for Cl albeit being the most versatile in that Cl is free 
from the well-known negative sign problem of QMC and the dimensional restriction of the 
tensor networkIn this view, extending the degrees of freedom accessible by Cl is a very 
important problem. 

In the held of quantum chemistry, the Cl wavefunction of an A^-electron system is ex¬ 
panded as 

M 

( 1 ) 


T (xi • • • xn) = ^ (u) 

in terms of the molecular orbitals (MOs) (x), which are represented as a linear combination 
of the atomic orbitals (AOs) (j)a, as 


M 


V’i (x) = ^ UiaCpa {x) , 


( 2 ) 


a=l 


using elements of an orthogonal matrix U. The Cl coefficients which are elements 

of an antisymmetric tensor of order A^jd] and dimension M, are varied to minimize the 
total energy and thus the full-CI calculation is, in practice, hampered by the degrees of 
freedom that grow combinatorially with N. Usually, the Cl series is truncated to make the 
computation tractable; however, when starting from the Hartree-Fock (HF) wavefunction by 
taking the permutation tensor in place of and the Hartree-Fock (HF) orbital 

in place of tjji, the Cl series usually converges very slowly |l|. The convergence does 
not, in general, speed up drastically even when starting from multiple Slater determinants. 
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Instead, use of localized MOs in place of the canonical HF orbitals is known to be effective for 


this purpose [^, and, recently, signihcant speed-up was achieved by using non-orthonormal 
Slater determinants comprised of non-orthonormal MOs that augment the HF orbitals jo]. 
It is therefore expected that Cl can be made more practicable by using optimal MOs for 
this purpose. 

Apart from Cl, the many-body perturbation theory such as the coupled cluster (CC)[^ 
has been developed to collect inhnite correction terms, thereby enabling the correction of the 
total energy of the water molecule, for example, within the error of 0.53 mHartree by taking 
up to triple excitations (CCSDT). However, the computational cost scales as 0{N^), which 
is still prohibitively high for application to many systems of interest. Thus, accelerating 
the convergence of the Cl series should open up the possibility to treat, with unprecedented 
accuracy, a whole new class of problems in electronic structure theory of strongly correlated 
systems. 

The key to the development that we present here is in a compact representation of the 
antisymmetric tensor With this in mind, we proposed in our previous workjs] to 

represent as a product of and a symmetric tensor and then to decompose 

the latter into a series of Kronecker products of a set of vectors such that 

K 

(3) 


r=l 


Such decomposition of an order-iV tensor into a linear combination of rank-1 tensors is known 
in the literature as the canonical decomposition { q] or the parallel factor decomposition 10]. 
In our approach, which we call the symmetric tensor decomposition Cl (STD-CI), variational 
parameters are the set of vectors c and the orthogonal matrix U and thus the degrees of 
freedom are KM + M {M — 1) /2. The number of order-1 tensors (i.e., the rank of the 
tensor decomposition iF ^), was found to be relatively few for small molecules (H 2 , He 2 , 
and LiH) and the Hubbard tetramer (FIG. 1-4). This suggests effectiveness of the tensor 
decomposition in compactly describing the wavefunction with a rapidly converging series 
(STD series) when the molecular orbitals are optimized together. Note that each term in 
the STD series contains all possible Slater determinants, which means that all the terms 
of the Cl series are regrouped differently to form different terms of the STD series. The 
computational cost thereby required is proportional to when using the algorithm 


to handle the permutation tensor developed in Ref. 


8j, so that the effectiveness of STD- 
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Cl depends on the the rank of the decomposition K required to accurately express the 
wavefunction. When applied to larger systems, however, we have found that the rank K 
needed for sufficient convergence is increased in general. In addition, the computation suffers 
from loss of signihcance in the floating-point arithmetic as will be shown below. Therefore, 
an improvement in the tensor decomposition is clearly desirable. 

In this context, one may use a more elaborate tensor decomposition, such as the tensor 
train network (TTN) suitable for one-dimensional svstemsllll.ll2| or the projected entangled 
pair states (PEPS) extended to two-dimensional svstemsjlSl]. Indeed, the former approach 
was applied to several molecular svstemsjl^. However, here we take a different approach: 
within the canonical decomposition technique, the degrees of freedom are extended by using 
different non-orthonormal MO coefficients for different terms in the STD series. That is, we 
use a general complex matrix for the MO coefficient although it was described as clUia in 
the original STD-CI scheme, extending thereby the degrees of freedom to KM^. We will call 
our new scheme the extended STD (ESTD). As will be shown below, each term in the STD 
series corresponds to a different Hartree-Fock-Bogoliubov (HFB) wavefunction 151]. Since the 
HFB wavefunction is called alternatively the antisymmetrized geminal power (AGP) 

ESTD may also be referred to as a linear combination of AGP. 

The AGP wavefunction has been used for variational calculations in a different context. 
In quantum chemistry, a single AGP (i.e., K = 1) has been used as the trial wavefunction 


to obtain the potential energy surface of small molecules 


(VMG) calculations, single AGP 22-25] or multi-AGP 26 


2l| . In variational Monte Garlo 


27j is multiplied by a correlation 


factor of the Jastrow type to form a Jastrow-AGP (JAGP) trial wavefunction; the opti¬ 
mization is done only for the Jastrow parameters with the MO coefficients kept unchanged 
from the initial HE orbitals or the initial natural orbitals constructed by performing a small 
preparatory GI calculation. In nuclear physics, multi-quasiparticle states are generated from 
a single AGP and are linearly combined to describe the wavefunction. This computational 
scheme is called the generator coordinate method (GGM) 28j]. The total-energy scheme 
of GGM was formulated by Onishi and Yoshida29|], and the mathematical structure was 


studied recently by other groups 


30 


3l| . Gontrary to those earlier works, we fully optimize 


the MO coefficients in our ESTD scheme without introducing the correlation factor. In the 
present scheme, we use an algorithm to decompose the permutation tensor into products of 
the second order permutation tensors allowing thereby to reduce the computational scaling 
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to K^M^. Contrary to STD-CI, this scheme does not require orthonormalization of AOs 
so that one can take advantage of the AO basis that is spatially localized and chemically 
comprehensible. We will show below that the STD series can be drastically shortened for 
the H 2 O molecule and Hubbard tetramers. This may open up the possibility to interpret 
strongly correlated systems in terms of MOs. In the next section, we introduce the outline of 
our formalism. We restrict our study to systems with even number of electrons throughout 
this paper. 


II. FORMALISM 


In the original STD-CI, we have applied the canonical decomposition to the symmetric 
part of the Cl coefficient (Eq. ([3])), so that the wavefunction Eq. ([T]) is rewritten using Eq. 
(EJ as 

M 

T (xi • • • Xtv) = ^ (a^l) • • • 0ajv {^ n ) (4) 

ai 

with the Cl coefficient becoming 


A 


(i\ 


M K 

E E4' 

»’=1 


’ ’ ’ ^ino-Ni 


(5) 


where K is the rank of the symmetric tensor decomposition and c[ is a set of complex 
vectors. We will take the convention for the subscript such that i, j,... correspond to MOs 
and a,b,... to AOs. In the extended STD-CI (ESTD), we expand the degrees of freedom 
by replacing the product with a general complex matrix such that 

KM K 


A, 


CLl---(lN 


E E 

r=l 


g. . m ...77^ 


E^ 

r=l 


ai •••aiv * 


( 6 ) 


Using the fact that the permutation tensor of order N for even N can be decomposed into 
a product of the second order permutation tensors as 

1 


(A/2)!2^/2 
1 


(cr) 




(tGSjv 




*3*4 '-iN-l'I-N J 1 


(7) 


( 8 ) 


(A/2)!2^/2 

where Sn is the symmetric group of degree N and A is the antisymmetrizer, Eq. ([6]) can be 
represented as 

(9) 


A 


a\ 


- AIM 7 ^' •••' 7 '’ ) 

(A/2)!2'^/2 ' la^ai /ajv-iajv* 
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using the antisymmetrized geminal 


ll, = Y^^nUlU],. 


( 10 ) 




This shows that our ESTD is an extension of the antisymmetrized geminal power (AGP) 
[is, 1^ in that AGPs are linearly combined. In the second-quantization form, the state 


vector is given by 


with 


K 


I't) = E 


r=l 


M 


17’') = 


(Ar/2)!2^/2 

t 


'^abCicl 


6=1 


N/2 


| 0 ) 


= exp 




t J 

a^b 


ab 


| 0 ) 


( 11 ) 

( 12 ) 

(13) 


iJV /2 


where the subscri^ means taking a coefficient of degree N/2 from the Hartree-Fock- 
Bogoliubov state 1^. This indicates a formal similarity of the ESTD wavefunction with 
that of the generator coordinate method (GGM) 28|], and thus one can take advantage of 


the formulas developec 


Onishi and Yoshida 


29 


(^2 l^^i) = exp 


for GGM. For example, we can use the matrix elements derived by 
such as 

Mr log (l (14) 

^ _ iN/2 

(16) 


= pf (1 

('''“IMIV) = (r^^^|^)^pf(i + 7-V‘t) 


{f^\c\clcdCc\-i^^) = ( 2 
+t 


1 


(16) 


iJV /2 


1 + 


7 


Ti 


cd - 


1 + 

1 


db 


7 


l- 2 t. 


ba 


pf (i 


, ( 17 ) 


tN/2 


1 + 

where 7 I corresponds to the Hermitian conjugate of 7 and pf denotes the Fredholm Pfaffian 
which is the square root of the Fredholm determinant 


pf (1 -|- = V det (1 -|- 7 '’ 2 t 7 ’'if). 


(18) 


It was shown in Ref. 321] that pf (l -|- 7 ’'^ 17 ’’R) is a polynomial of degree M/2 since roots 
of the characteristic polynomial of a product of two antisymmtric matrices are pairwise 
degenerated. 
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The total energy can be obtained from the matrix element of the Hamiltonian ("H) 

|7^1) = ^ hab {Y^\clca I 7 ’'') + ( 7 ''' I cMcdCc b’’'), (19) 

ab abed 

where h corresponds to the sum of the kinetic energy and the external potential and V 


is the Coulomb repulsion. Contrary to standard GCM schemes, we explicitly ob 


coefficient of in Eqs. (ll41ll71) instead of using the particle projection method(29| 


similar projection method 


29|, 


ain the 


32|. A 


32| can be applied to the HFB state to adapt the wavefunction 


to the spatial and spin symmetries, which is necessary to obtain accurate ground state 
wavefunctions, but we do not use it in the present study to compare ESTD and full-CI 
without the symmetry adaptation. We also explicitly obtained the derivatives of the energy 
with respect to 7 , with which to optimize the parameters U[^ using the Broyden-Fletcher- 
Goldfarb-Shanno algorithm (BEGS) SSHSTj. We emphasize that ESTD allows the handling 
of non-orthogonal AOs by replacing 7 ^^ by where s is the overlap matrix of 

AOs, although orthogonality is assumed for simplicity in the present formulation. 

When using general complex matrices for in practical numerical works, the computa¬ 
tion is sometimes impossible to carry on because of the loss of significance in the floating¬ 
point arithmetic. To avoid this problem, we make the geminal complex unitary as well as 
antisymmetric so that all the eigenvalues are equal to unity in the absolute values although 
the convergence speed is slightly affected. Note that the geminal can be made unitary by 
the following redefinition of noticing that the second order permutation tensor e is real 


antisymmetric, we can apply the canonical tranformation 
0 A \ / 0 A ^/2 \ / Q 


e = U^ 


U = U^ 


38| by which 
I \ I 0 A72 


U, 


( 20 ) 


-AO/ \ A72 0 } \ -I 0 I \ A72 0 

where A and I are, respectively, real diagonal and identity matrix, G is a real orthogonal 
matrix, and the superscript T means taking the matrix transpose: then the geminal (Eq. 
flTOj) 1 becomes 

Y = (u'-f ( ^ 

where has been redefined as 

0 A72 

A72 0 


Ur 


( 21 ) 




UU^ 


( 22 ) 
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Since t/'’ was introduced in Eq.(|6]) as a general matrix to be determined variationally, this 
redefinition does not affect the resulting formulas. In this way the geminal becomes complex 
unitary when restricting to be complex unitary. We call such restricted ESTD the 
antisymmetrized unitary geminal powers (AUGP) hereafter. The variation of the unitary 
matrix can be done using the Cayley transform 

U = {1 + X)-\l - X) (23) 

with a skew-Hermitian matrix X. 

III. RESULTS 

We have applied ESTD and AUGP to the water molecule and half-filled Hubbard 
tetramers and compared the total energy with the STD-CI and the full-CI calculations. 
In the water molecule, we set the geometry condition to 

O = (0, 0, 0), i/ = (-1.809, 0, 0), (0.453549,1.751221, 0) (24) 

in atomic units. The calculation was done using the basis set STO-3G. The minimal basis 
set is sufficient for comparison between the different schemes, although larger basis sets will 
be required for highly quantitative calculations. The Hubbard tetramer has a tetrahedral 
geometry, and all the sites are connected by the transfer integral of unity. The Hubbard U 
was taken as 1, 10^, and 10“^. 

Figures [H |2] and [3] show the total energy referred to the full-CI calculation, or the residual 
error vs. the tensor rank of the decomposition. The convergence is much faster for ESTD 
than for STD-CI. For larger U cases, the residual error is initially almost insensitive to the 
rank and the value amounts to 1-10 when rescaled in unit of U. Then the residual error 
drops suddenly when K is larger than 2 for [/ = 10^ and 3 for U = 10^. This might indicate 
that the anti-ferromagnetic (AF) ground state cannot be described even at a qualitative 
level when the number of parameters is unreasonably small. For U = 10“^, the total energy 
is almost the same among ESTD, STD-CI, and full-CI when the rank is 6. We postulate 
that the parameter set for STD-CI happens to become suitable for describing the AF limit 
only at 77 > 6. 

Table m shows detailed comparison of the total energy calculated by AUGP and full-CI 
for the Hubbard tetramer with U = 10^. The calculation was done using double precision 



K (number of terms) 


FIG. 1: Energy error of 4 sites U =1.0 Hubbard model for ESTD (solid line with circle) and STD-CI 
(broken line with circle). 



K (number of terms) 


FIG. 2: Energy error of 4 sites 17=100.0 Hubbard model for ESTD (solid line with circle) and 
STD-GI (broken line with circle). 

with K = 3 . The residual error is 2.2 x 10“^^ while the error was around 10“^^ for ESTD 
(Fig. ([2])); in either scheme the error is close to the double precision limit. 

Figure m shows the residual error of the total energy of the H 2 O molecule. ESTD shows 
strikingly faster convergence than STD-CI. In the ESTD and STD-CI calculations, we used 
quadruple precision since serious loss of significance occurred with double precision. Because 
of this problem, we could not easily obtain accurate results using K larger than 4 for ESTD, 
and thus we could not reduce the residual error below 10“^ Hartree. Figure [5] and Table UTI 
show the result obtained by using AUGP with double precision. In obtaining the residual 
error, we used the total energy from full-CI calculation as the reference, but the result 
is almost the same when using the complete active space self-consistent held (CASSCF) 


9 








FIG. 3: Energy error of 4 sites C/=10000.0 Hubbard model for ESTD (solid line with circle) and 
STD-CI (broken line with circle). 


Method 

Total energy 

AUGP 

Exact 

-0.119880248924950 

-0.119880248946222 


TABLE I: Total energy of the Hubbard tetramer with U = 100 obtained by the AUGP with K = 
3. The result obtained by full-GI is also shown. 


calculation performed using GAUSSIAN09 package j39|. The figure shows that the residual 


energy is reduced almost linearly in log-scale with respect to K, indicating an exponential 
convergence. The table shows that the error is 6.4 x 10“^ Hartree with K = 4 and is reduced 
to 7.2 X 10“^ Hartree with K = 10. Although AUGP yields slower convergence than ESTD, 
the AUGP is very stable. It is not difficult to obtain the optimal parameters LTj even when 
the calculation is started from random initial values. 


IV. CONCLUSION 

To efficiently represent the GI wavefunction, we have applied the canonical decomposition 
algorithm to the symmetric part of the GI coefficient. Therein the molecular orbitals are 
fully optimized, without imposing the orthonormality condition, differently for different 
terms in the decomposition series. The computational scheme, which we call the extended 
symmetric tensor decomposition (ESTD), is equivalent to the linear combination of the 
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FIG. 4: Residual error of the total energy of the H 2 O molecule obtained by ESTD (solid line with 
circle) and STD-CI (broken line with circle). 


Method 

Total energy 

AUGP, K = 1 

-61.508355 

AUGP, K = 2 

-72.985876 

AUGP, A = 3 

-73.457483 

AUGP, A = 4 

-75.006008 

AUGP, A = 6 

-75.011870 

AUGP, A = 8 

-75.012385975 

AUGP, A = 10 

-75.012425100 

Exact (ours) 

-75.012425818 

Exact (Gaussian) 

-75.012425839 


TABLE II: Total energy (in unit of Hartree) of H 2 O obtained by AUGP. Eor comparison, the full- 
GI calculation was done using our own code and the GASSGE calculation using the GaussianOD 
package. 

Hartree-Fock-Bogoliubov states or the linear combination of the antisymmetrized geminal 
power (AGP). By this, we can rearrange the full-CI series into the canonical decomposition 
series (STD series). By applying ESTD to the water molecule and Hubbard tetramers, we 
found that the total energy rapidly converges well within ten terms {K = 10) for the STD 
series. The ESTD calculation was found to be numerically unstable when the number of 
electrons is increased as large as 10 because of the loss of signihcance in the floating-point 
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FIG. 5: Residual error of the total energy of H 2 O obtained by AUGP 


arithmetic. This problem could be avoided by restricting the MO coefficients to be complex 
unitary although the convergence speed with respect to tensor rank K was slightly affected. 
The computational cost of ESTD scales as where M is the number of MOs. The 

result suggests that the AGP-based scheme is a promising computational tool for quantum 
systems. In this context, it will be an important target of future study to clarify how K 
depends on the complexity and the size of the system. Our calculation was done without 
parallelization, but an acceleration by a factor of can be expected because of the almost 
independent nature of the computation. Further acceleration is expected by applying the 
tensor decomposition scheme to the two-electron integrals as done in Ref. 40| ; such technique 
may possibly reduce the scaling to 
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